library(GEOquery)
library(tidyverse)
library(tidySummarizedExperiment)
library(tidybulk)
library(ggpubr)

gset <- getGEO("GSE83452", GSEMatrix =TRUE, getGPL = TRUE, AnnotGPL = TRUE)

gset <- gset[[1]]

probe_gene <- gset@featureData@data |>
  filter(GB_ACC != '') |>
  as_tibble()

final_geneid <- probe_gene$GB_ACC |>
  clusterProfiler::bitr(fromType = 'ACCNUM', toType = 'SYMBOL',OrgDb = 'org.Hs.eg.db') |>
  filter(str_detect(SYMBOL,'GSDM|DFNA5|ICERE')) |>
  left_join(probe_gene, join_by(ACCNUM == GB_ACC)) |>
  select(2:3)

smmex <- makeSummarizedExperimentFromExpressionSet(gset)
smmex

max(smmex@assays@data$exprs)

smmex <- smmex |>
  select(1:3, characteristics_ch1.1, characteristics_ch1.2) |>
  mutate(disease = str_remove(characteristics_ch1.1, 'liver status: '),
         intervention = str_remove(characteristics_ch1.2, 'type of intervention: '))

smmex_coarse <- smmex |>
  filter(.feature %in% final_geneid$ID)

smmex_coarse |>
  as_tibble() |>
  filter(intervention == 'Diet' & disease != 'undefined') |>
  left_join(final_geneid, join_by(.feature == ID)) |>
  mutate(disease = fct_relevel(disease, 'no NASH')) |>
  ggplot(aes(disease, exprs, color = disease)) +
  geom_jitter(height = 0, width = .1) +
  stat_summary(geom = 'errorbar', fun.data = 'mean_sd', width = .3) +
  stat_summary(geom = 'crossbar', fun = 'mean', width = .5) +
  stat_compare_means(method = 't.test', comparisons = list(c('NASH','no NASH'))) +
  expand_limits(y = 0) +
  facet_wrap(~SYMBOL, scales = 'free') +
  scale_y_continuous(expand = expansion(mult = c(.1,.1))) +
  scale_color_manual(values = c('blue','red')) +
  theme_pubr()

smmex_coarse |>
  as_tibble() |>
  left_join(final_geneid, join_by(.feature == ID)) |>
  select(2,3,7,8) |>
  write_csv('NASH_GSDMX_GSE83452.csv')

# GSE130991 ------
gset <- getGEO("GSE130991", GSEMatrix =TRUE, getGPL = TRUE, AnnotGPL = TRUE)

gset <- gset[[1]]

probe_gene <- gset@featureData@data |>
  as_tibble() |>
  filter(seqname != '---') 

probe_accession <- probe_gene |>
  separate(gene_assignment, into = c('acc','symbol','other'), sep = ' // ') |>
  filter(symbol != '---') |>
  select(1,symbol)

final_geneid <- probe_accession |>
  filter(str_detect(symbol, 'GSDM|DFNA5|ICERE'))

smmex <- makeSummarizedExperimentFromExpressionSet(gset)

smmex_meta <- smmex@colData |>
  as_tibble() |>
  select(2, characteristics_ch1.2,
         characteristics_ch1.5,
         characteristics_ch1.6,
         characteristics_ch1.12:characteristics_ch1.15)


smmex_meta |> map(some, \(x)str_detect(x, 'diabetes')) 

notreat <- smmex_meta |>
  mutate(.keep = 'unused',
         statin = str_remove(characteristics_ch1.2, 'statin: '),
         bmi = str_remove(characteristics_ch1.5, 'bmi: ') |> as.numeric(),
         bmi2 = str_remove(characteristics_ch1.6, 'bmi: ') |> as.numeric(),
         bmi = case_when(is.na(bmi) ~ bmi2, .default = bmi)) |>
  filter(!(str_detect(characteristics_ch1.12, 'diabetes')|
             str_detect(characteristics_ch1.13, 'diabetes')|
             str_detect(characteristics_ch1.14, 'diabetes')|
             str_detect(characteristics_ch1.15, 'diabetes')))

notreat |>
  count(statin)

no_statin <- notreat |>
  select(geo_accession, statin, bmi) |>
  filter(statin == 'No Statin')

max(smmex@assays@data$exprs)

expr_mat <- smmex@assays@data$exprs |>
  as_tibble(rownames = 'ID') |>
  right_join(final_geneid)

expr_coarse <- expr_mat |>
  pivot_longer(where(is.numeric)) |>
  right_join(no_statin, join_by(name == geo_accession)) |>
  mutate(group = if_else(bmi > 30, 'BMI>30', 'BMI<=30'))

expr_coarse <- expr_coarse |>
  mutate(symbol = ifelse(symbol == 'DFNA5', 'GSDME', symbol))

expr_coarse |>
  ggplot(aes(group, value, color = group)) +
  geom_jitter(height = 0, width = .2) +
  stat_summary(geom = 'errorbar', fun.data = 'mean_sd', width = .3, color = 'black') +
  stat_summary(geom = 'crossbar', fun = 'mean', width = .5, linewidth = .3, color = 'black') +
  stat_compare_means(method = 't.test', comparisons = list(c('BMI>30','BMI<=30'))) +
  #expand_limits(y = 0) +
  facet_wrap(~symbol, scales = 'free') +
  scale_y_continuous(expand = expansion(mult = c(.1,.1))) +
  scale_color_manual(values = c('blue','red'), label = c('BMI<=30 (n=26)', 'BMI>30(n=549)')) +
  theme_pubr() +
  labs(title = 'GSE130991')

expr_coarse |>
  select(-ID) |>
  write_csv('Statin-diabete_GSDMX_GSE130991.csv')

# mouse data ----------
## GSE69306 ============
gset <- getGEO("GSE69306", GSEMatrix =TRUE, getGPL = TRUE, AnnotGPL = TRUE)

gset <- gset[[1]]

probe_gene <- gset@featureData@data |>
  as_tibble() |>
  filter(str_detect(Description,"asdermin") & !str_detect(Description,"pseudo|like") | ORF == 54722) |>
  mutate(suffix = str_to_lower(str_extract(Description, '(?<=\\s).+')),
         symbol = ifelse(ID == '54722_at', 'Gsdme', str_glue('Gsdm{suffix}'))) |>
  select(ID, symbol) 

smmex <- makeSummarizedExperimentFromExpressionSet(gset)

smmex_meta <- smmex@colData |>
  as_tibble() |>
  select(2, source_name_ch1) |>
  separate(source_name_ch1, into = c('group','tissue'), sep = '_') |>
  filter(str_detect(group, 'iet'))

smmex_meta

expr_mat <- smmex@assays@data$exprs |>
  as_tibble(rownames = 'ID') |>
  right_join(probe_gene) |>
  pivot_longer(where(is.numeric), names_to = 'geo_accession') |>
  right_join(smmex_meta) |>
  mutate(group = fct_relevel(group, 'Normal Diet'))

expr_mat |> count(tissue)

expr_mat |>
  filter(tissue == 'JEJUNUM') |>
  ggplot(aes(group, value, color = group)) +
  geom_jitter(height = 0, width = .2) +
  stat_summary(geom = 'errorbar', fun.data = 'mean_sd', width = .3, color = 'black') +
  stat_summary(geom = 'crossbar', fun = 'mean', width = .5, linewidth = .3, color = 'black') +
  stat_compare_means(method = 't.test', comparisons = list(c('High-fat Diet','Normal Diet'))) +
  #expand_limits(y = 0) +
  facet_wrap(~symbol, scales = 'free') +
  scale_y_continuous(expand = expansion(mult = c(.1,.2))) +
  scale_color_manual(values = c('blue','red')) +
  theme_pubr() +
  labs_pubr() +
  labs(subtitle = 'GSE69306', title = 'RNA expression in HFD mice jejunum',
       y = 'expression level')
  
expr_mat |>
  select(-ID) |>
  write_csv('DigestiveTract_GsdmX_GSE69306.csv')

## GSE63175 =======
gset <- getGEO("GSE63175", GSEMatrix =TRUE, getGPL = TRUE, AnnotGPL = TRUE)

gset <- gset[[1]]

probe_gene <- gset@featureData@data |>
  as_tibble()

View(probe_gene)

probe_gene <- gset@featureData@data |>
  as_tibble() |>
  filter(str_detect(GENE_SYMBOL,"Gsdm[a-d]|Dfna5") & GENE_SYMBOL != 'Gsdmcl1') |>
  mutate(GENE_SYMBOL = ifelse(GENE_SYMBOL == 'Dfna5', 'Gsdme', GENE_SYMBOL)) |>
  select(ID, GENE_SYMBOL)

smmex <- makeSummarizedExperimentFromExpressionSet(gset)

smmex_meta <- smmex@colData |>
  as_tibble()

smmex_meta <- smmex@colData |>
  as_tibble() |>
  separate(source_name_ch1, into = c('group','tissue'), sep = '_') |>
  filter(str_detect(group, 'LF|HFC')) |>
  select(2, group, characteristics_ch1.3)

smmex_meta

smmex_meta <- smmex_meta |>
  mutate(day = str_extract(characteristics_ch1.3, '\\d+'),
         time_unit = if_else(str_detect(characteristics_ch1.3, 'Day'), 1, 7),
         day = as.numeric(day) * time_unit) |>
  select(geo_accession, group, day)
  
expr_mat <- smmex@assays@data$exprs |>
  as_tibble(rownames = 'ID') |>
  right_join(probe_gene) |>
  pivot_longer(where(is.numeric), names_to = 'geo_accession') |>
  right_join(smmex_meta) |>
  mutate(group = if_else(group == 'LF', 'LFC', group),
         group = fct_relevel(group, 'LFC'))

expr_mat |>
  summarise(expr = list(value), .by = c(GENE_SYMBOL, group, day)) |>
  pivot_wider(names_from = group, values_from = expr) |>
  rowwise() |>
  mutate(pval = t.test(HFC, LFC)$p.value) |>
  filter(pval < .05)

expr_mat |>
  ggplot(aes(day, value, color = group)) +
  stat_summary(fun = 'mean', geom = 'path') +
  stat_summary(geom = 'errorbar', fun.data = 'mean_sd', width = .3) +
  #stat_compare_means(method = 't.test', comparisons = list(c('HFC','LFC'))) +
  #expand_limits(y = 0) +
  facet_wrap(~GENE_SYMBOL, scales = 'free') +
  scale_y_continuous(expand = expansion(mult = c(.1,.2))) +
  scale_color_manual(values = c('blue','red')) +
  theme_pubr() +
  labs_pubr() +
  labs(subtitle = 'GSE63175', title = 'RNA expression in HFD mice liver',
       y = 'expression level')

expr_mat |>
  filter(GENE_SYMBOL == 'Gsdmd') |>
  ggplot(aes(day, value, color = group)) +
  stat_summary(fun = 'mean', geom = 'path', linewidth = 2) +
  stat_summary(geom = 'errorbar', fun.data = 'mean_se', width = 1) +
  #stat_compare_means(method = 't.test', comparisons = list(c('HFC','LFC'))) +
  #expand_limits(y = 0) +
  scale_y_continuous(expand = expansion(mult = c(.1,.2))) +
  scale_color_manual(values = c('blue','red')) +
  theme_pubr() +
  labs_pubr() +
  geom_text(aes(x = 63, y = 3.1, label = '*'), inherit.aes = FALSE, size = 10) +
  geom_text(aes(x = 105, y = 3.1, label = '*'), inherit.aes = FALSE, size = 10) +
  labs(subtitle = 'GSE63175', title = 'Gsdmd RNA expression in HFD mice liver',
       y = 'expression level')

expr_mat |>
  select(-ID) |>
  write_csv('liverTime_GsdmX_GSE63175.csv')

### submitter's logFC result ------------
read_delim('GSE63175_fold_change_Liver_Day1_to_W18.txt.gz')


## GSE30534 =====
gset <- getGEO("GSE30534", filename = 'GSE30534_series_matrix.txt.gz', GSEMatrix =TRUE, getGPL = TRUE, AnnotGPL = TRUE)

probe_gene <- gset@featureData@data |>
  as_tibble()

View(probe_gene)

probe_gene <- gset@featureData@data |>
  as_tibble() |>
  filter(str_detect(GeneSymbol,"Gsdm[a-e]|Dfna5") & GB_ACC != '') |>
  select(ID, GeneSymbol) |>
  mutate(GeneSymbol = ifelse(GeneSymbol == 'Dfna5h', 'Gsdme', GeneSymbol))

smmex <- makeSummarizedExperimentFromExpressionSet(gset)

smmex_meta <- smmex@colData |>
  as_tibble()

View(smmex_meta)

smmex_meta <- smmex_meta |>
  mutate(day = str_extract(title, '\\d+') |> as.numeric(),
         tissue = str_remove(characteristics_ch1.2, 'tissue: ')) |>
  select(geo_accession, day, tissue)

smmex_meta

expr_mat <- smmex@assays@data$exprs |>
  as_tibble(rownames = 'ID') |>
  right_join(probe_gene) |>
  pivot_longer(where(is.numeric), names_to = 'geo_accession') |>
  right_join(smmex_meta)

expr_mat$value |> max()

expr_mat |>
  ggplot(aes(day, value, color = tissue)) +
  stat_summary(fun = 'mean', geom = 'path') +
  stat_summary(geom = 'errorbar', fun.data = 'mean_sd', width = .3) +
  #stat_compare_means(method = 't.test', comparisons = list(c('HFC','LFC'))) +
  #expand_limits(y = 0) +
  facet_wrap(~GeneSymbol, scales = 'free') +
  scale_y_continuous(expand = expansion(mult = c(.1,.1))) +
  scale_color_manual(values = c('blue','red','green','orange','#007FFF','magenta')) +
  theme_pubr() +
  labs_pubr() +
  labs(subtitle = 'GSE30534', title = 'RNA expression in HFD mice\'s tissues',
       y = 'expression level')

expr_mat |>
  select(-ID) |>
  write_csv('TissueTime_GsdmX_GSE30534.csv')

## GSE216442 ========
gset <- getGEO("GSE216442", GSEMatrix =TRUE, getGPL = TRUE, AnnotGPL = TRUE)

gset <- gset[[1]]

probe_gene <- gset@featureData@data |>
  as_tibble()

View(probe_gene)

probe_gene <- gset@featureData@data |>
  as_tibble() |>
  filter(seqname != '---') 

probe_accession <- probe_gene |>
  separate(mrna_assignment, into = c('acc','db','symbol'), sep = ' // ') |>
  mutate(GeneSymbol = str_extract(symbol, '\\(.+\\)')) |>
  filter(str_detect(GeneSymbol,"Gsdm[a-e][^l]|Dfna5")) |>
  select(ID, GeneSymbol) |>
  mutate(GeneSymbol = ifelse(str_detect(GeneSymbol, 'Dfna5'),
                             'Gsdme',
                             str_extract(GeneSymbol, '\\w+') ))

smmex <- makeSummarizedExperimentFromExpressionSet(gset)

smmex_meta <- smmex@colData |>
  as_tibble()

View(smmex_meta)

smmex_meta <- smmex_meta |>
  mutate(tissue = str_remove(tissue.ch1, ' homogenate'),
         gender = gender.ch1, diet = diet.ch1) |>
  filter(diet != '45% HFD') |>
  select(geo_accession, gender, tissue, diet)

smmex_meta

expr_mat <- smmex@assays@data$exprs |>
  as_tibble(rownames = 'ID') |>
  right_join(probe_accession) |>
  pivot_longer(where(is.numeric), names_to = 'geo_accession') |>
  right_join(smmex_meta) |>
  mutate(diet = fct_relevel(diet, 'control'))

expr_mat$value |> max()

expr_mat |>
  summarise(expr = list(value), .by = c(GENE_SYMBOL, group, day)) |>
  pivot_wider(names_from = group, values_from = expr) |>
  rowwise() |>
  mutate(pval = t.test(HFC, LFC)$p.value) |>
  filter(pval < .05)

expr_mat |>
  filter(tissue == 'cerebellar' & gender == 'female') |>
  ggplot(aes(diet, value, color = diet)) +
  geom_jitter(height = 0, width = .2) +
  stat_summary(geom = 'errorbar', fun.data = 'mean_sd', width = .3, color = 'black') +
  stat_summary(geom = 'crossbar', fun = 'mean', width = .5, linewidth = .3, color = 'black') +
  stat_compare_means(method = 't.test', comparisons = list(c('60% HFD','control'))) +
  #expand_limits(y = 0) +
  facet_wrap(~GeneSymbol, scales = 'free') +
  scale_y_continuous(expand = expansion(mult = c(.1,.2))) +
  scale_color_manual(values = c('blue','red')) +
  theme_pubr() +
  labs_pubr() +
  labs(subtitle = 'GSE216442', title = 'RNA expression in HFD female mice cerebellar',
       y = 'expression level')

expr_mat |>
  select(-ID) |>
  write_csv('anteCor&crbllr_GsdmX_GSE216442.csv')

## GSE63174 ======
gset <- getGEO("GSE63174", GSEMatrix =TRUE, getGPL = TRUE, AnnotGPL = TRUE)

gset <- gset[[1]]

probe_gene <- gset@featureData@data |>
  as_tibble()

View(probe_gene)

probe_gene <- gset@featureData@data |>
  as_tibble() |>
  filter(str_detect(GENE_SYMBOL, "Gsdm[a-e]|Dfna5")) |>
  filter(GENE_SYMBOL != 'Gsdmcl1') |>
  mutate(GENE_SYMBOL = ifelse(str_detect(GENE_SYMBOL, 'Dfna5'),
                             'Gsdme',
                             GENE_SYMBOL)) |>
  select(ID, GENE_SYMBOL)
  

smmex <- makeSummarizedExperimentFromExpressionSet(gset)

smmex_meta <- smmex@colData |>
  as_tibble()

View(smmex_meta)

smmex_meta <- smmex_meta |>
  mutate(group = fed.with.ch1,
         time = str_extract(time.point.ch1, '\\d+') |> as.numeric(),
         day = ifelse(str_detect(time.point.ch1, 'Week'), time * 7, time)) |>
  select(geo_accession, group, day) |>
  filter(!str_detect(group, 'KAL'))

expr_mat <- smmex@assays@data$exprs |>
  as_tibble(rownames = 'ID') |>
  right_join(probe_gene) |>
  pivot_longer(where(is.numeric), names_to = 'geo_accession') |>
  right_join(smmex_meta) |>
  mutate(group = fct_relevel(group, 'Low fat diet'))

expr_mat$value |> max()

expr_mat |>
  filter(day != 1) |>
  summarise(expr = list(value), .by = c(GENE_SYMBOL, group, day)) |>
  pivot_wider(names_from = group, values_from = expr) |>
  rowwise() |>
  mutate(pval = t.test(`Low fat diet`, `High fat diet`)$p.value) |>
  filter(pval < .05)

expr_mat |>
  ggplot(aes(day, value, color = group)) +
  stat_summary(fun = 'mean', geom = 'path') +
  stat_summary(geom = 'errorbar', fun.data = 'mean_sd', width = .3) +
  #stat_compare_means(method = 't.test', comparisons = list(c('HFC','LFC'))) +
  #expand_limits(y = 0) +
  facet_wrap(~GENE_SYMBOL, scales = 'free') +
  scale_y_continuous(expand = expansion(mult = c(.1,.2))) +
  scale_color_manual(values = c('blue','red')) +
  theme_pubr() +
  labs_pubr() +
  labs(subtitle = 'GSE63174', title = 'RNA expression in HFD mice hippocampus',
       y = 'expression level')

### separate genes
expr_mat |>
  filter(GENE_SYMBOL == 'Gsdmd') |>
  ggplot(aes(day, value, color = group)) +
  stat_summary(fun = 'mean', geom = 'path', linewidth = 2) +
  stat_summary(geom = 'errorbar', fun.data = 'mean_sd', width = 1) +
  geom_text(aes(x = 14, y = 3.1, label = '***'), inherit.aes = FALSE, size = 10) +
  #stat_compare_means(method = 't.test', comparisons = list(c('HFC','LFC'))) +
  scale_y_continuous(expand = expansion(mult = c(.1,.2))) +
  scale_color_manual(values = c('blue','red')) +
  theme_pubr() +
  labs_pubr() +
  labs(subtitle = 'GSE63174', title = 'Gsdmd RNA expression in HFD mice hippocampus',
       y = 'expression level')

expr_mat |>
  filter(GENE_SYMBOL == 'Gsdme') |>
  ggplot(aes(day, value, color = group)) +
  stat_summary(fun = 'mean', geom = 'path', linewidth = 2) +
  stat_summary(geom = 'errorbar', fun.data = 'mean_sd', width = 1) +
  geom_text(aes(x = 14, y = 3.1, label = '*'), inherit.aes = FALSE, size = 10) +
  geom_text(aes(x = 42, y = 3.1, label = '*'), inherit.aes = FALSE, size = 10) +
  #stat_compare_means(method = 't.test', comparisons = list(c('HFC','LFC'))) +
  scale_y_continuous(expand = expansion(mult = c(.1,.2))) +
  scale_color_manual(values = c('blue','red')) +
  theme_pubr() +
  labs_pubr() +
  labs(subtitle = 'GSE63174', title = 'Gsdme RNA expression in HFD mice hippocampus',
       y = 'expression level')

expr_mat |>
  select(-ID) |>
  write_csv('hippocampus_GsdmX_GSE63174.csv')
